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We present a semi-analytical model to simulate bidirectional reflectance distribution function (BRDF) 
spectra of a rough slab layer containing impurities. This model has been optimized for fast computation 
in order to analyze hyperspectral data. We designed it for planetary surfaces ices studies but it could be 
used for other purposes. It estimates the bidirectional reflectance of a rough slab of material containing 
inclusions, overlaying an optically thick media (semi-infinite media or stratified media, for instance gran¬ 
ular material). The inclusions are supposed to be close to spherical, and of any type of other material 
than the ice matrix. It can be any type of other ice, mineral or even bubbles, defined by their optical 
constants. We suppose a low roughness and we consider the geometrical optics conditions. This model is 
thus applicable for inclusions larger than the considered wavelength. The scattering on the inclusions is 
assumed to be isotropic. This model has a fast computation implementation and thus is suitable for high 
resolution hyperspectral data analysis. © 2015 Optical Society of America 

OCIS codes: (010.5620) Radiative transfer; (240.6490) Spectroscopy, surface; (110.4234) Multispectral and hyperspectral imaging 
; (240.5770) Roughness ; (290.7050) Turbid media ; (080.2468) First-order optics 
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1. introduction 

Hyperspectral imaging has become a major component in plane¬ 
tary surface observation since the past decades. Earth and other 
Solar system bodies are now observed in various spectral ranges 
at various resolutions and from various heights. 

As a photon come across a surface, it interacts in two major 
ways. It can be either absorbed or deviated (scattering, diffrac¬ 
tion, refraction). The objective of this radiative transfer model is 
to describe the interactions using a realistic surface description. 
In such descriptions, the reflectance of a surface is the result 
of multiple interactions, with multiple irregular interfaces of 
different materials. The exact resolution of the radiative transfer 
equations turns out to be a highly difficult and time consuming 
problem. This problem has been solved under certain hypoth¬ 
esis : if the characteristic size of the particle is much smaller 
than the wavelength, or if it is much bigger. In this study, we 
consider the geometrical optics domain, where the particles are 
much bigger than the wavelength. For example, in the visible 
and near infrared range (400 nm — 5 pm), we suppose that the 
average particle size does not fall below 10 pm. This is in general 
valid for planetary surfaces [1]. Ray tracing algorithms [2-5] that 
simulate the very complex paths of millions of photons through 
these surface can show very accurate results, but they depend on 
a huge number of parameters and are highly time consuming, se¬ 


riously limiting the the interpretation of extensive hyperspectral 
images. We aim at a radiative transfer model that is fast enough 
to be able to deal with a vast amount of data, such as planetary 
spectro-imaging databases. It is then necessary to make further 
simplifying assumptions that enable the formulation of approxi¬ 
mate analytic or semi-analytic solutions to the radiative transfer 
problem. A possible simplification is to consider that the radia¬ 
tive properties inside a media can be described statistically only 
using local mean properties of scattering and absorption [6-8]. 
The media is assumed to be homogenous at a mesoscopic scale. 
Another classical simplifying assumption considered in such 
problem is the two stream approximation [6, 7, 9]. It has been 
shown that under certain conditions, it did not affect too much 
the solution compared to more accurate studies, but simplifies 
greatly the calculations [10,11]. To describe the reflectance of a 
surface, one also have to consider the geometry of illumination 
and observation. In our approach, these photometric effects 
are modeled by the properties of the interface between the me¬ 
dia and the exterior. These properties of roughness can also be 
statistically described, using only one or a few parameters. Y. 
Shkuratov [8] and B. Hapke [7] developed analytical radiative 
transfer models for granular media, that are able to simulate the 
bidirectional reflectance of various granular surfaces. 

If the media cannot be described as homogenous, it is possible 
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to consider it pieciwise continuus, constituted of homogeneous 
strata. It is the case for example in the atmospheres, or stratified 
surfaces. A family of models describe the radiative transfer in 
stratified media, such as the DISORT algorithm [12]. In these 
discrete-ordinate modelisations, each layer is considered ho¬ 
mogenous, and the total reflectance is calculated iteratively, by 
adding the contribution of each layer. This method has also been 
adapted to the study of the ocean-atmosphere coupled system 
with a rough surface [13]. 

Starting from Hapke model, improving it, and combining it 
with a multi-layer method [12], S. Doute has also developed a 
model for stratified granular surfaces [14]. Using the same strat¬ 
egy, we developed a semi-analytical radiative transfer model for 
a compact layer (solid matrix containing inclusions) overlaying 
an optically thick granular layer. This two layers approach does 
not require an iterative DISORT-like method, but only adding 
coupling formulas. It is founded on three major assumptions 
: (i) the geometric optics conditions are observed, (ii) the me¬ 
dia is piecewise continuous and (iii) the inclusions are close to 
spherical and homogeneously distributed in the matrix. 

Model overview 

We decompose the reflectance into two distinct contributions 
: specular and diffuse. We chose the Hapke [15] probability 
density function of orientations, as it well describes the statistic 
distribution of slopes in the approximation of small angles. We 
consider a collimated incident radiation, at an incident angle i. 
We estimate the specular contribution, considering the geometry 
and the surface description. The specular reflection of rough 
surfaces have been studied in various cases [15-19]. We use 
the same general idea of these methods, describing the rough 
surface as constituted of multiple unresolved facets. The specu¬ 
lar contribution will result from the integration of the specular 
reflections on the facets, in the solid angles considered the 
light source and the detector) as described in figure la. 

Then we estimate the diffuse contribution. The total reflec¬ 
tion coefficient at the first rough interface, that determines the 
amount of energy transmitted to the slab, is obtained by inte¬ 
grating specular contributions in every emergent direction, at a 
given incidence. We consider that the first transit through the 
slab is anisotropic (collimated), and that there is an isotropisa- 
tion at the second rough interface ( i.e. when the radiation reach 
the semi-infinite substrate). For the refraction and the internal 
reflection, every following transit is considered isotropic. The 
diffuse contribution is obtained using an analytical estimation 
of Fresnel coefficients [14, 20], and a simple statistical approach. 
The contribution of the semi infinite substrate is estimated using 
Hapke model [21]. Finally, we consider that the slab is under a 
collimated radiation from the light source, and under a diffuse 
radiation from the granular substrate. We compute the resulting 
total bidirectional reflectance using adding doubling formulas 
(figure lb). 

2. SURFACE ROUGHNESS - FACETS DISTRIBUTION 

The first step is to describe the roughness of the surface. We 
consider that it is composed of N facets that are not resolved, 
with N 1. These facets' orientations follow a probability den¬ 
sity a ($, £), where & is the zenital angle between the normal to 
the facet and the local vertical direction, and £ is the azimutal 
angle. To make our approach as general as possible, we chose to 
describe the surface as randomly rough. Such a roughness has 
already been widely studied (see for example [15-19,22] and the 
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Fig-1. Illustration of the radiative transfer in a rough slab. 

(a) Radiative transfer for a slab ice layer only. Anisotropic 
transit are represented in red, and named with a prime. On the 
top left: illustration of the reflections and transmission at the 
first interface, used in the calculations of variables S' e and fact, 
0'. (b) Illustration of the adding coupling. The granular and 
slab layers are artificially separated in this figure to help the 
understanding of the coupling. 
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reference cited in these papers). These studies show that a slope 
distribution, with tan 8 that is close to Gaussian is a good de¬ 
scription of the surface. Such a description combines simplicity 
and efficiency reproducing the photometric variations. For the 
sake of simplicity and because it is widely used in the literature, 
we chose the following probability distribution function [15] : 


a (*,£) = 


n L tarn 


= exp 


tan 2 8 
n tan 2 8 


sec 2 8 sin 8 (1) 


where 


relations, taking into account the geometry variations within 
one pixel and the statistics of the slopes a (8, £). 

A. Specular conditions for one facet 

For one facet Af to satisfy the specular reflection conditions, 
then its normal Wf must respect (see figure 2): 

E = -I + 2{I\W f )-W f (3) 

where the operator (|) represents the scalar product. If we ex¬ 
press Eq. 3 in the ( x,y,z ) coordinates,: 


2 r" 

tan0=— / 2 a ($) tantfd# (2) 

n Jo 

It is supposed that the azimutal distribution is uniform. The 
angle 8 representing the mean slope angle completely charac¬ 
terizes the facets' orientations and the surface roughness. This 
slopes distribution considers the cases of small 8. Practically, 
the threshold of validity can be determined depending on the 
level of tolerance (see Figures 5 and 6. The expression of a (8, f) 
could be adapted in the future to extend the study to any type 
of terrain, as discussed in section A. 



Fig. 2. The local coordinate system ( x,y,z) is centered at the 
slab surface, with the z axis vertically upward and the x axis 
horizontally toward the sun. I and E are respectively the inci¬ 
dent and the emergent directions. Wf is the normal to the facet 
A f . 

3. SPECULAR REFLECTANCE 

In this part we determine the specular contribution to the re¬ 
flectance at the detector. We first establish the relations between 
the orientation ( 8 spec , Cspec) of a facet that is in specular condi¬ 
tions and the geometry of observation (z, e, ip), where i is the in¬ 
cidence angle, e the emergence and ip the azimuth (see Figure 2). 
The total specular contribution is obtained by integrating these 


! x e = 2 (sin i sin 8 spec cos L, sp e C + cos 8 spec cos i) sin 8 spec cos £ spec - sin i 
ye — 2 (sin z sin 8 spec cos £s pec T cos 8 spec cos z) sin 8 spec sin ^spec 
z e = 2 (sin i sin 8 spec cos £ spec + cos 8 spec cos z) cos 8 sp ec — cos i 

(4) 

where x e = sin e cos ip, y e = sin e sin ip and z e = cos e are the 
coordinates of the emergent vector E in the (x, y, z) frame. This 
leads to : 


tan 2 tispec 
COS E,spec 
Sin £spec 


_ sin 2 z’+sin 2 e+2 sin z sin e cos ip 

(cosz+cose) 2 

= sHTt^v (l + tan 2 8 S p ec ) — cosz) 

sin e sin ip 
(cos z+cos e ) tan 

(5) 


B. Expression of the specular reflectance for one pixel 

We consider a pixel of area A formed of N facets of same area 
Af, orientated according the probability density a (8, f), as de¬ 
tailed in section 2, with N»l. The number of facets satisfy¬ 
ing the specular reflection conditions defined in section A will 
be JJ nc N a ( 8 spec , l , spec ) d (8, £), where H c is the set of values 
{8, Q satisfying 5 within the range of observation geometries. 
Indeed, there is a range of different geometry of observation 
within one instruments' pixel. Let Xc be the emergence varia¬ 
tions within a pixel. The facet orientation satisfying the specular 
conditions for this geometry is ( 8 spec , £spec) given by Eq. 5. At 
the incidence z, there is a range of orientations that satisfy specu¬ 
lar condition within a pixel that is centered at (8 spec ,^ spec ) and 
has the size 8 (8, £). <5 (8, £) is determined using a function 

g. : (S. ('■♦) = « 

V'/y \g2( e y) = c 

that transforms (e,tp) into ($,£), for the incident angle i. Not 
every facet that satisfy these conditions will send energy to 
the captor. Indeed, the roughness of the surface introduces a 
shadowing of the scene : some facets will not receive incident 
light, or will not be visible by the captor, or both. A shadowing 
factor S / must be introduced at this point. Let N spec be the 
number of facets that satisfy both the geometrical condition 
defined in section A and the visibility condition : 

N S p ec — If N a ($ S pec, Zspec)S' {i,e,xp,8) d{8,Z) (6) 

J •j'Hc 

with S ; (z, e, ip,8) a shadowing factor that depends on the geom¬ 
etry of observation and the roughness of the surface [15]. Each 
one of these N S pec facets receives an incident power P; and send 
back a reflected power P r : 

Pi = FA f cos (7) 
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F being the incident power flux ( e.g. : Solar flux) in the radi¬ 
ation direction, Aj cos ^ y J the projection of the facet in the 

plane orthogonal to the incident radiation, and Tf ^ 

Fresnel reflection coefficient in energy at the phase angle a 1 , 
Tf = R 2 ± (a) + Rj (a). As a' does not depend on the facet ori¬ 
entations, all these specular reflections will result in a specular 
power Pspec = NspecPr, thus : 



angle of the source. Then the total specular contribution within 
one pixel will be : 

Rspec ( i , e, ip) = R'spec {h e , V) | d et Jg e (i, ip ) | sin i d i dip 

(13) 

where det Jg e (z, ip) I is the Jacobian of the function 


* 


FT ('-«-*) 

W \S2 V 7 ) =CJ 


that transforms ( i,xp ) into (#,£), for the given emergence angle 
e. 


Pspec (i,e,ip,9) = ff NFA 
J Jl-Lc 


(Aspect Zspec)S' d(#,£) (9) 


The reflectance factor R is the ratio between the bidirectional 
reflectance r of the surface and the one r L of a perfectly lam- 
bertian surface , thus R = 7r^-j. The bidirectional reflectance 
r is the ratio between the radiance L of the surface and the 
collimated incident power, perpendicularly to the incident di¬ 
rection. Thus r = y, with L = QAcose 7 ^ being the illuminated 
surface and Cl c the solid angle subtended by a pixel. Finally, 


R. 


spec 


= 71- 


— 7T 


‘-•spec 

F cos i 


= 71 


1 spec 

Cl c AF cos i cose 


, thus : 


r r NAf COS 

Rspec (i,e,y,e) = jj Hc n ncAcosicos S' {i,e,ip,e) 


($spec/ £spec ) d(tf,o do) 


where O c is the solid angle subtended by an instrument's pixel. 
A is the sum of the horizontal projection of all the facets : A = 
NAf (cos#). The term (cos#) is included in the shadowing 
function S (i, e, ip, 9) described by B. Hapke [15]. Thus we can 
simplify Eq. 10 as 


R. 


spec 


!L 


-S ( i,e,ip,9 ) 


H c Oc cos i cos e 

T f ^" 2 "^ a (#s pec/Cspec) d(#,£) (11) 


d(#,£) is derived from the integration angles e and ip that are 
the emergence and azimuth angles. There is a bijection between 
Oc and He because 5 admits a unique solution for every (e, ip). 
Considering that the incidence angle i is a constant, we can 
rigorously express Rl pec as : 


R 


i 

spec 


{i,e,ip) 




cos i cos e 


S (i, e, i p,8) Tf 



a(#s pec^spec) det 7g. (e, i/’) dedip 


( 12 ) 


| det /g e (e, ip) | is the Jacobian of the function g ; . This expres¬ 
sion 12 assumes that the incidence i is a constant. In reality, the 
light source is almost never completely collimated, but ranges 
inside a solid angle (e.g. : the Solar disk). Let Q s be the solid 


4. DIFFUSE REFLECTANCE 

We consider a two layers model, with a slab overlaying a semi¬ 
infinite granular substrate. The collimated radiation from the 
sun is transmitted to the slab with an external reflection coeffi¬ 
cient S' e (the prime here represent the anisotropy). We suppose 
an isotropisation at the second interface. The slab is modeled 
as a compact isotropic and homogeneous matrix. It contains 
inclusions that are close to spherical and not identical to the ma¬ 
trix. The inclusions are the main contributors to the scattering 
of radiation in the layer. They are distributed homogeneously 
in the matrix. The determination of the Fresnel coefficients at 
the interface matrix/inclusion or inclusion/matrix is a key to 
estimate the transmission and reflection factors of the layer. An 
internal and external reflection coefficient S,* and S^ for each 
type of inclusion k must be defined. 

In this part we describe the radiative transfer in the media. 
First we will characterize the transmission of light into the slab. 
By energy conservation this is equivalent to calculating the to¬ 
tal reflected power which, normalized by the incident energy, 
stands for the reflexion coefficient (section A). Then we will de¬ 
scribe the scattering of light by the inclusion during the transfer 
through the slab. This requires the calculation of the external 
and internal reflection coefficients of these inclusions (section B). 
Once the basic optical properties of the inclusions are known, 
we can consider fluxes of energy within the whole slab that will 
be governed by the radiative properties of the slab (section E). 
Solving this radiative transfer problem within the slab with an 
upper and lower optical interfaces will give the overall reflec¬ 
tion and transmission factors of the slab (section F). Finally the 
radiative interaction of the two layers (substrate and slab) are 
considered and solved by adding doubling leading to the final 
result (section G). 


A. Reflection coefficients for the slab 

Anisotropic case Let S' be the external reflection coefficient in 
a collimated case (interface atmosphere/ice matrix). It corre¬ 
sponds to the ratio between the incident power P, and the total 
reflected power, in every direction P,! of . The total reflected power 
can be estimated integrating the specular contributions for every 
emerging direction, at the given incidence angle i. 


Si 


JJln dPspec 
AF cos i 


(14) 


d P S p ec being the specular contribution in a given geometry. Us¬ 
ing Eq. 9, the expression of S' e becomes 



(#spec/ £,spec) S (i, C, ijl, 9 j 

cos i 


d (#,0 

(15) 
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where H is the set of values taken by & and £ throughout the 
integration. Exactly like in section B, d($, £) is derived from the 
integration angles e and ip that are the emergence and azimuth 
angles. There is now a bisection between B and T~L, B being the 
superior hemisphere that is the domain of variation of e and 
ip. Considering that the incidence angle i is a constant, we can 
express S' e as : 


A r 2n COS 1 

(i) 

|r/| 


1 a ($ spec/ £ spec') S (/, 6,Xp,0^ 

Jo J 0 




COS / 


det /g. (e, ip) de dtp (16) 


where |det J s (e,ip)\ is the Jacobian of the function 


Finally, the external reflection coefficient S ek at the interface 
matrix/inclusion is : 

7T 

2 

s ek = J [r 2 ^ (a) +R| (a)] e^ a "P(i-cos a ) cosn;sinada (20) 
o 

This differential absorption effect is already taken into account 
in the expression of R± (ot) and R|| (a) in the case of the internal 
reflection at the interface inclusion/matrix. 

C. Fresnel coefficients 

The Fresnel reflectivities for perpendicular and parallel polar¬ 
ization with respects to the propagation plane, for an incidence 
angle ot, R± (a) and R|| (ot), are derived from Snell's law (see [21] 
Hapke, 2012, sect. 4.3, pp.46-60): 


g : M _ (*■ <'•*>-' 

ypj \S2( e >f) = C 

that transforms (e, tp) into ($, f), for the incident angle i. 

The internal reflection coefficient S' in a collimated case at 
the interface ice matrix/atmosphere is not considered as we 
suppose an isotropisation of the radiation at the second interface 
(ice/granular regolith). 

Isotropic case In the isotropic case, the internal reflection co¬ 
efficient S, is obtained integrating the Fresnel equations at the 
surface for the all geometries: 

7T 

2 

Si = j ^R“y ( a ) + («)J cos a da (17) 

0 


R,l (a) = 


^( a)= (co sa-^+0| 

(cos ot + Q\ ) 2 + Q\ 

[( n 2 — fc 2 ) cosa — Gi] 2 + [2nfccosa — Q 2} 2 
[( n 2 — fc 2 ) cos a + Qi] 2 + [2rcfccosa + Q 2] 2 


( 21 ) 


( 22 ) 


using n = n i n 2 +hkz and ^ _ tiih nzh , m + ifci and «2 + ik 2 
the complex refractive indexes of the media considered. 

1 


Qi = 


c 2 — 1 

e 2 - 2 


| n 2 — fc 2 — sin 2 aj + (ir — fc 2 — sin 2 a j + 4 n 2 k' 


2 . .1 

1 "I 

2 

) +4n 2 fc 2 



(23) 


■ 

' / „ \ 2 ^ 

1 -1 

2 

— ^n 2 — fc 2 — sin 2 aj + 

n 2 — fc 2 — sin 2 ct j + 4n 2 fc 2 



(24) 


where R± (a) and (a.) are Fresnel reflectivites for perpendic¬ 
ular and parallel polarization with respects to the propagation 
plan, for an incidence angle ot and will be detailed later. 

The external reflection coefficient S e is estimated the same 
way : 

7T 

2 

S e = J (a) + R 2 (a)j cos a da (18) 

0 


D. Integration 

S. Chandrasekhar showed (see Chandrasekhar, 1960, sect. 22, 
pp. 61-69 [20]) that many radiative transfer integration can be 
approximated using the Gauss quadrature formulae. If / (p) is 
a polynomial of order 2m — 1, then 

[ f (R) = E c if Ui) (25) 

-1 i =1 


B. Reflection coefficients for the inclusions 

In the case of a spherical inclusion of the type fc, the internal 
reflection coefficient Sj k is obtained in the usual fashion inte¬ 
grating the Fresnel equations (see [21] Hapke, 2012, sect. 5.4.4, 
pp.78-95) 

7r 
2 

S, 7 t = J \r?i (tx.) + R| (a)j cos a sin a da (19) 

0 

For the estimation of the external reflection coefficient S ek , 
a differential absorption factor is taken into account. Indeed, 
as we deal with inclusions in an absorbing matrix, the parallel 
rays we consider in the integration touch the inclusion after 
different optic paths. For a ray that touches the inclusion with 
an incidence ot, the differential path length in the matrix is v = 
p k cos ot, where p/ c is the radius of the spherical inclusion. Thus 
the differential absorption factor is e ““P'/ 1 cosa ) / where a m is 
the absorption coefficient of the matrix. Writing the matrix's 
optical index n m + i k m , the dispersion relation gives a m = 


where are the zeros of the Legendre polynomials 

Pi ,..., P m of order 1,..., m, and c\,cj are the associated Christof- 
fel numbers : 



Equation 25 is exact if / (/() is a polynomial of order 2m — 1. 
When f (p) is not a polynomial, then the quadrature formulae 
gives an approximation that converges to the exact value when 
m —> 00 . The order m of the approximation directly governs 
its quality. We estimate analytically the internal and external 
reflection coefficients in the isotropic case S, an S e using the 
roots of the 32^* order Legendre's polynomial and the associ¬ 
ated Christoffel's numbers as detailed in [23]. We use a simple 
change of variable to transform the integration interval from 
[0, j] into [—1,1]. All the integrations are performed using the 
Gauss quadrature formulae, except S' e and R S p ec . In these cases, 
the integration being a double one, we cannot use the Gauss 
quadrature. We chose after numerical tests an adaptive grid and 
the rectangle method. 
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E. Radiative properties of a slab containing inclusions 

We suppose an homogeneous distribution of isotropic inclu¬ 
sions inside the slab. The inclusions type is noted by k, with 
Nj different types, defined by different geometrical and optical 
properties. 


E.l. Proportions of inclusions 

We define the slab compactness y c as the volume of the ice 
matrix per unit of volume. We also define Af the total number of 
inclusions per unit of volume, and A 4 the number of inclusions 
of the type k per unit of volume. The proportion of each type of 
inclusion is P k = 4^. Immediate geometrical calculations give : 


Af = 


3 (7c - 1) 

^LkU p kPk 


(27) 


E.2. Cross sections 

We suppose close to spherical inclusions. The scattering effi¬ 
ciency for a sphere has been described by B. Hapke ([21] Hapke, 
2012, sect. 5.6, pp.95-99, eq 5.52a) in his equivalent slab model. 
For an inclusion of the type k : 

Qsk = S ek + (l-S ek ) ^J s S Q ik ®ik (28) 

where S[ k and S ek are respectively the internal and external re¬ 
flection coefficients of an inclusion expressed in equations 19 
and 20, and ©,j. is the internal transmission coefficient of and 
inclusion. In the two stream approximation, and assuming the 
isotropy of the phase function of the internal scatterers in an 
inclusion ([21] Hapke, 2012, sect. 6.5, pp.122-144, Eq. 6.26) the 
expression of <dj k can be reduced simply to : 

r ik + exp o k3 /a lk ( a ik +s ik )) 

®ik = -7- . \ (29) 

1 + rj k exp [ p k \/ a lk ( a + s ; / c ) j 

aj k being a type k inclusion's absorption coefficient, Sj k its scat¬ 
tering coefficient, and 


rik = 


1 - 
1 + 




a ik 

a ik J r^ik 

a ik 

a ik+ s ik 


The scattering cross section cr sk for one inclusion is 


(30) 


a travel of length di/, the probability p\ for a photon to meet an 
inclusion and be scattered is: 


P\= 1 - exp (-Af (c t s ) ^2^ 


(33) 


The probability p 2 that this photon has not been absorbed by the 
matrix before is: 

p 2 = exp (~a m dv) (34) 

Thus the probability p s for a photon to be only scattered per unit 
of length is: 


Ps = ^ ex P (~ a mdv) 


1 — exp ( —Af {cr s ) ^ ^ dv 


that becomes for an infinitesimal length dv: 

p s =Jf(cr s ) ln7f +o(l) 
7c -1 


(35) 


(36) 


equally, the probability p 3 for a photon to be absorbed or scat¬ 
tered by an inclusion throughout dv is: 


p 3 = exp (—fl m dv) 


1 — exp ( —Af {a e ) J * 11 —j-dv 


(37) 


and the probability to be absorbed by the matrix during dv is: 

Pi = 1 — exp (—a,„dv) (38) 

so the probability of extinction p e per unit of length is: 


Pe = 


dv 


1 - exp ( - (Jf (cr e ) 


+ a m dv 


when di' is close to 0 , it becomes: 


Pe = Jf ( cr e ) 


In 7c 
7c -1 


In 7c 
7c -1 


+ a m + o ( 1 ) 


(39) 


(40) 


Finally we obtain the single scattering albedo of a slab containing 
inclusions dividing p s by p e : 


Af {cr s } 


* <^> + lh77« 


(41) 


^sk — &kQsk ( 31 ) 

where cr* is the geometrical cross section : a \ = np Let (a s ) be 
the mean cross section of the inclusions : 

Ni 

(Vs) = L P k cr sk (32) 

k=l 

We do the approximation of geometric optics, so the extinction 
cross section a ek correspond to the geometrical cross section a k . 

E.3. Single scattering albedo and optical thickness 

The single scattering albedo u> of an absorbing and scattering ob¬ 
ject is defined as the ratio of the total amount of power scattered 
to the total amount of power removed to the wave (absorbed or 
scattered). We propose a simple statistical approach to express 
the single scattering albedo of a unit of volume of slab contain¬ 
ing inclusions. We use the same method as [21] Hapke, 2012, sect. 
7.4, pp.158-169, but we modify the medium description. After 


Equation 39 gives the expression of the optical depth t of a slab 
with inclusion: 

r = (cr e ) ^ Lj + a,^j v (42) 

F. Diffuse reflectance and transmission factors of the contam¬ 
inated slab 

F. 1. Diffuse reflectance of a slab under collimated illumination 

In this section, we suppose that the slab is under a collimated 
radiation. As in section 2, we suppose that the surface is con¬ 
stituted of N unresolved facets that have a slope distribution 
given by the probability density function a (#,£). Each facet 
will receive an illumination at an incidence i f depending on its 
orientation. We consider in that case that the first transit in the 
slab is collimated and will transmit rays of light into the slab at 
different inclinations. Our goal at this point is to determine the 
mean transmission path D' trough a slab of a given roughness 
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9 and a given thickness D. For a facet with the orientation de¬ 
fined by (9, £) using Snell-Descartes law gives D' = fact($, £)D, 
with: 


| ~ h cos *'+ cos& (h cos */ - \A “ S ( X " cos2 if )) | 

(43) 

i | being the incidence angle on the facet (i.e. the angle be¬ 
tween the facet's normal and the incident radiation). Using 
basic trigonometric relations gives [21]: 

cos if = sin i sin 9 cos £ + cos 9 cos i (44) 

We consider that only the first transit in the slab is anisotropic. 
The internal absorption factor for the first anisotropic transit 
0' will depend on the mean length Lf of this transit, with 
D' =factD , and 

2n f 

fact = I j fact (#,£) a (#,£) d#d£ (45) 

o o 

thus 


F.2. Diffuse reflectance of a slab under isotropic illumination 

In this section we suppose that the slab is under an isotropic 
radiation. Indeed, at the lower interface, it is illuminated isotrop¬ 
ically from below by the substratum. Rg and Tg have their usual 
expressions in this case : 


Rq — Se + 


(1 - Se) S/0 2 (1 - Si) 
1 - (0S,-) 2 


(52) 


= 0(1-S e )(l-S,-) 
1 - (0S,) 2 


(53) 


G. Bidirectional reflectance of a contaminated slab overlaying 
a semi-infinite granular media 

In realistic conditions, a slab will receive a collimated radiation 
from the solar disk, and a diffuse radiation from the granular 
medium underneath. There is a coupling between the two layers, 
illustrated on Figure 1. Using adding doubling formulas [14], 
we can express the total diffuse reflectance of the slab over a 
granular substrate as : 


fact = 


n 2 tan 6 


2n 2 

IJ 


_ tan z & 

e n tan^ o sec 2 i? sin $ 


- sb cos f + cos * (it cos b - v /l “^( 1 “ cos2; /)) 


d9df 

(46) 


and 

r m + exp (■ -D' y a m (a m + p s ) ] 

0' = -y- ' (47) 

1 +r m exp y-D' sja m (a,„ + p s )J 

where p s is given by Eq. 36, r m = hi , and co is given by 

Eq. 41. 

The internal absorption factor for an isotropic transit is ([21], 
Hapke, 2012 Eq. 6.26) 


0 = 


r m + exp ^-2 Dsja m (a m +s)J 
1 + r m exp r—2D sja m (a m + s) j 


(48) 


R D iff = K + T o T ors E ( R o r s) n 


n =0 


— 1^0 + 


'Tq'ToD 
1 - Ro r s 


(54) 


where r s 


1 - 

1 + 


is the lambertian reflectance of the sub¬ 


strate [14]. oj s is the single scattering albedo of the granular 
substrate. The last step is to simulate the diffuse contribution 
for one measurement. The total reflectance (BRDF) of the sur¬ 
face measured by the instrument is the sum of the specular and 
diffuse contributions : 


Rtot — Rspec + RDiff (55) 

where R S pec is determined by Eq. 13 and RDiff by Eq. 54. 

5. DISCUSSION 


Every following transit is considered isotropic. As illustrated 
on figure 1, we can express the reflectance of the slab under a 
collimated radiation Rq as 


R'' = S' + (1-S')0'S / ©(1-S / ) 


i + E ( 0S ;) 


n =1 


(49) 


S, being the internal reflection coefficient of the slab. The term S' e 
represents the integration over the sky of the specular reflectance, 
and the other represents the diffuse reflectance. Thus we can 
express the diffuse reflectance of the slab as 


K = 


(1 - S') &Sj 0(1-Sf) 
1 - (©Si) 2 


(50) 


The diffuse transmission of the slab under a collimated radiation 
Tq is obtained the same way : 


v = 0'(1-S') (1 — Si) 

0 i - (©s ,) 2 


(51) 


A. Energy conservation 

At the first interface We checked the conservation of the energy 
at different points in the model. We first checked it at the first 
interface, as it contains a complex numerical integration. To test 
the conservation of the energy at the first interface, we force the 
value of the Fresnel's reflection coefficient r f in Eq. 16 to one. 
Thus, all the energy is supposed to be sent back, and we have to 
obtain Q = 1 to have the energy conserved, where 


A r 2n cos 

Q= / 

-'0 * / 0 


a(# s ,£s)S(i,e, ip, 9) 
cos i 


det7g. (e,i p) 


de dip 

(56) 


that is Eq. 16 where the Fresnel reflection coefficient is put to 
one. Figure 3 shows the value Q as a function of the incidence 
angle. Different roughness parameters 9 were tested, ranging 
from 9 = 0.01° to 9 = 45°. Only values ranging from 9 = 0.15° 
to 6 = 3.5° are displayed on Figure 3. This test illustrates the 
dependance of the validity of the model on both the incidence 
angle and the roughness parameter. 
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Energy conservation versus incidence 



Incidence angle (°) 


Fig. 3. Q as a function on the incidence angle i, when forcing 
the value of the Fresnel's reflection coefficient r f in Eq. 16 to 
one. A value of one means that the conservation of energy is 
respected. This figure shows that in the cases of a roughness 
parameter ranging from 8 = 0.15° to 8 = 3.5°, the energy is 
fairly well conserved for incidences below 85°, and a rough¬ 
ness parameter below 8 = 2.5°. Thus, this model will not be 
applicable to very high incidences values, and roughness over 
8 = 2.5° 



Energy conservation versus incidence 
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Fig. 4. Q as a function on the incidence angle i, where £,■ is 
the energy sent back from the surface, integrated over the 
hemisphere, and £,■ the incident energy. 


For the complete model To test the conservation of the energy 
for the whole model, we first have to set the complex value of 
the optical constant of the slab and the granular substrate to 0, to 
make the surface non absorbent. Then we integrate the energy 
sent back toward the sky. This energy must equal the incoming 
energy. To test this practically in the model, we set the sensor's 
angular aperture to a value that is equal to the integration step. 

Figure 4 show the value of Q = ^ J'q f 0 2 Rtot cose sinede dip. 
Practically, the energy is conserved if this integral equals 1. In¬ 
deed, the energy conservation gives 


/ 


LA cos e 
FA cos i 


dQ 


1 


(57) 


where £ is the surface radiance (W.m _2 .sr _1 ) A is the surface of 
a pixel, f is the incident flux in the incident direction (W.m~ 2 ) 
and i and e are the incidence and emergence angles. The relation 
R — 7iYcosi between the reflectance factor and the radiance 
brings 

1 2 n f 

— I J R to t cos e sin e de dip = 1 (58) 

o o 

The symmetry of the model in azimuth leads to the quantity Q 
displayed in Figure 4. This figure shows that it is mostly the 
roughness parameter 8 and the incidence angle i that control the 
validity of the model. 

Figure 5 shows the error in the energy conservation in percent, 
as a function of the roughness parameter 8 and the incidence i. 
This gives the range of validity of the model according to a given 
tolerance. Roughness parameters larger than 8 = 11°, always 
exhibiting error larger than 10 %, are not represented. For small 
slab real optical index (i.e. close to one), these errors decrease. 


Error in the conservation of energy 



Fig. 5. Error in the energy conservation as a function of 6 and 
i. 
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A. 1. Slope distribution 

As mentioned in section 2, this model is limited to the case of 
small 6. Figure 5 gives a quantification of that limitation. This 
is mostly due to the fact that the probability density function 
a (9, f) that defines the repartition of slopes only makes sense if 
././(i? g a (#, £) d# d£ = 1, which means that 


1 

2n 



= 1 


(59) 


— / \ 

or that the value of I = f " —exp ( — tan2 f- ) sec 2 d sin i9 di? 
must be equal to 1. Figure 6 shows the value of I versus 
the roughness parameter 9. The function a (#,£) only makes 
sense as a probability function if I = 1. For values of rough¬ 
ness larger than 9 = 2°, I begins to fall to values below one. 
In a further development, we could extend the applicability 
of the model by defining a new probability density function 
ANorm (#,£), that would be the normalization of the function 



Fig. 6. I versus 9. For values of 9 larger than 9 = 2°, the 
probability density function for the slopes begins to drop. At 
6 = 10°, I = 0.957. 


B. Behavior of the model 
B.l. Specular reflection 

This model is built to simulate observations. Thus, the specular 
spot characteristics will depend not only on the illumination 
divergence and the geometry, but also on the observation de¬ 
vice. This model is designed to be adaptable to both conditions. 
Figure 7 shows a zoom on the specular spot for a water ice slab 
at 1 pm with a roughness parameter 9 = 0.5°, illuminated at an 
incidence angle i = 50° illuminated with different light sources, 
and observed with two distinct detectors. In the first case (in Fig¬ 
ure 7a), the surface is illuminated with a light source that has an 
aperture of 0.4° and observed with a sensor that has an aperture 
of 4.2°. It represents the conditions of a laboratory measurement 
with the instrument described in [24]. In the second case (in Fig¬ 
ure 7b), the surface is illuminated with a light source that has an 
aperture of 0.2° and observed with a sensor that has an aperture 
of 6.92.1CP 20 . It represents the conditions of a measure with the 
OMEGA imaging spectrometer instrument orbiting the planet 
Mars [25]. Both cases represent actual measurement situations. 


As shown on Figure 7, both the amplitude and the shape of the 
specular spot depend on the characteristics of the illumination. 

B.2. Influence of the parameters 

To give a feeling on how the model behaves according to the 
different parameters, we chose a set of parameters, and plotted 
the dependence of the reflectance to the variation of one parame¬ 
ter around this first set. We chose arbitrary optical constants for 
the matrix and inclusions. We chose for the matrix n = 1.3 and 
k = 1.1CP 3 , that are approximately the values for water ice at 
270 K, and at the 1 pm wavelength. We selected as our standard 
set of parameter a 10 mm thick slab layer containing 1000 ppmv 
of 100 pm wide inclusions, overlaying a semi infinite granular 
layer of the same nature that the matrix. We tested the behavior 
of the model for various types of inclusions. We describe two 
types of behavior. The first type is when the absorption coef¬ 
ficient of the inclusions is smaller than the one of the matrix. 
This includes the particular case of a matrix contaminated with 
bubbles. It is illustrated in the figure 8,10 and 9 with the green 
curves. In this case, the real part n of the optical index of the 
inclusions has very little influence. The second case is when the 
absorption in the inclusion is bigger than in the matrix (blue and 
red curves). In this case, the real optical index of the inclusions 
plays an important role. 

Figure 8 shows the dependence of the reflectance on the thick¬ 
ness of the slab layer in different cases. In the case of an uncon¬ 
taminated slab, the reflectance approaches 0 when the thickness 
increases. On the contrary, the reflectance of a slab containing 
inclusions will saturate at a value depending on the properties 
of the impurities. In the first case of a low absorption in the 
inclusions (green curve) the reflectance is higher than the re¬ 
flectance of an uncontaminated layer whatever thickness the 
slab has. On the other case of high absorption in the inclusions, 
when the slab layer is thin, then the reflectance is lower than 
one of an uncontaminated slab, but as thickness increases, the 
value saturates, due to the scattering of light by the inclusions. 
This value gives an idea of the penetration depth of the light 
into a contaminated slab layer (i.e. the depth from which the 
layer becomes optically thick). Figure 9 show the dependence 
of the reflectance on the radius of inclusions in the slab layer. 
This illustrates the scattering properties of the inclusions. The 
reflectance factor of an uncontaminated at this geometry is ap¬ 
proximately R = 0.372. In this figure, the volumetric proportion 
of inclusions in the matrix is held constant, so as the grain size 
increases, the number of inclusions per unit of volume decreases. 
Thus the scattering power decreases as well as a function of the 
grain size. In the case of inclusions with higher absorption than 
the matrix, at some point, the grain size reaches a value where 
the absorption in the inclusions becomes more efficient than the 
scattering effect, and the reflectance falls below the reference 
value of the uncontaminated slab (around 10 pm for the blue 
curve and 20 pm for the red one). Then, the decreasing proba¬ 
bility of encountering an inclusion when grain sizes become too 
high make the reflectance approach the value of the uncontami¬ 
nated slab. Indeed, when the grain size of the inclusion equals 
the thickness of the layer (at the extreme right of the plot), the 
probability of encountering one, knowing the the volumetric 
proportion is 1000 ppmv becomes very low and the influence of 
the inclusions become negligible. 

Figure 10 shows the evolution of the reflectance in different 
cases where scattering or absorption dominates. On Figure 10a, 
the absorption is the dominant effect. In the case of inclusions 
with a higher absorption coefficient than the matrix, this makes 
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Fig. 7. Zoom on the specular spot for a water ice slab at 1 pm 
with a roughness parameter 6 — 0.5°, (a) illuminated at an 
incidence angle i = 50° with a light source that has an aperture 
of 0.4°, and observed with a detector that has an aperture of 
4.2°, such as in the conditions of a laboratory measurement 
with the instrument described in [24], and (b) illuminated at 
an incidence angle i = 50° with a light source that has an aper¬ 
ture of 0.2°, and observed with a detector that has an aperture 
of 6.92.10~ 2 °, such as in the conditions of a measure with the 
high resolution spectro imaging imaging instrument OMEGA 
orbiting Mars described in [25]. 


i = 50.0° e= 10.0“ ij/= 180.“ p= lOQG.ppm g= lOO^tm 



Fig. 8. Reflectance factor of a slab of water ice containing 
various types of inclusions, at a wavelength A = 1 pm, as a 
function of the thickness of the slab layer, other parameters 
fixed. Black curve is the uncontaminated reference, and col¬ 
ored curves represent different optical indexes of inclusions. 



1.0 10.0 100.0 1000.0 10000.0 
Radius of the inclusions (urn) 


Fig. 9. Reflectance factor of a slab of water ice containing 
various types of inclusions, at a wavelength A = 1 pm, as a 
function of the grain size of the inclusions, other parameters 
fixed. Colored curves represent different optical indexes of 
inclusions. When the absorption coefficient is higher in the 
inclusions than in the matrix (blue and red curves), there is a 
competition between scattering and absorption. The dominant 
effect depend on the grain size. The reflectance factor of an 
uncontaminated at this geometry is approximately R = 0.372. 
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i = 50.0“ e= 10.0" rp= 180." h= lO.mm g= lOO.jum 




i = 50.0° e= 10.0“ rp= 180“ h= 10.mm g= lO.jum 




Fig. 10. Reflectance factor of a slab of water ice containing 
various types of inclusions, at a wavelength A = 1 pm, as 
a function of the volumetric proportion of inclusions, other 
parameters fixed. Colored curves represent different optical 
indexes of inclusions, (a) absorption is the dominant effect 
for blue and red curves, (b) scattering is the dominant effect 
for red curve, and scattering and absorption are of the same 
order of magnitude (blue curve). The reflectance factor of an 
uncontaminated at this geometry is approximately R = 0.372. 


i = 50.0“ e= 10.0“ ip= 180.° 



Fig. 11. Reflectance factor of a 20 mm thick slab of water ice 
containing 1000 ppm of 100 pm wide inclusions as a function 
of the roughness parameter 9, other parameters fixed. The 
optical indexes of the inclusions are n = 1.1 and k = 1.10~ 3 . 
The reflectance decrease as the roughness increase. This is 
due to the fact that a bigger roughness means more facets in 
specular conditions and thus less energy inserted into the 
system. The roughness parameter have a smaller impact on 
the diffuse reflectance than the other parameters of the model. 


the reflectance drop when the proportion of inclusions increases 
(blue and red curves). For the green curve, both absorption and 
scattering contribute in increasing the reflectance of the slab, 
thus it increases with the proportion of inclusions. Figure 10b 
shows more complexity. It is the same as Figure 10a except that 
the grain size of the inclusions is 10 pm instead of 100 pm. The 
green curve still represents the case of a lower absorption in the 
inclusion. Absorption and diffusion both tend to increase the 
reflectance, thus it increases with the proportion of inclusions. 
The red curve represents a case of higher absorption in the in¬ 
clusions, when the scattering contribution is dominant. In this 
case, as diffusion limits the penetration of light into the layer, the 
reflectance increases with the proportion of impurities. The blue 
curve represents the limit case where diffusion and absorption 
contributions are of the same order of magnitude, leading to 
strong non-linear behavior. 

Figure 11 shows the dependence of the reflectance of a 20 mm 
thick slab of water ice containing 1000 ppm of 100 pm wide inclu¬ 
sions on the roughness parameter 9. The optical indexes of the 
inclusions are n = 1.1 and k = 1.10~ 9 . The diffuse reflectance of 
a slab decreases as the roughness increases. A bigger roughness 
means a bigger diversity in the slope distribution at the surface. 
This leads to an increased number of facets satisfying the specu¬ 
lar reflection conditions defined in section A. Finally, less energy 
is inserted into the surface, and the diffuse reflectance is smaller. 
The dependance of the diffuse reflectance on the roughness if 
smaller compared to the dependance on the other parameters. 
This can be attributed to the relatively small range of values 
of 9. On the contrary the roughness parameter 9 has a strong 
influence on the specular contribution. 
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6. CONCLUSIONS 

We developed a radiative transfer model to simulate the bidi¬ 
rectional reflectance of a rough slab with inclusions. Typical 
calculation time is 1.10 _2 s per spectrum, considering 10000 
wavelengths. Nevertheless, it can vary greatly depending on the 
sets of parameters desired. 

Most of the constituting elements of this model have already 
been numerically [14] or experimentally validated [21] but we 
adapted them to build a new model. In this study we tested 
numerically the conservation of energy and characterized the 
domain of validity of the model. We conducted sensibility stud¬ 
ies in the case of a matrix containing only one type of inclusion. 
This shows the complexity and non linearity of the model with 
respect to its parameters. The sensibility study in the case of 
several types of inclusions were not conducted because of the 
great number of parameters. The experimental validation will 
be conducted in a following paper. 

This model is designed to analyze massive hyperspectral 
data in the planetary science domain. In our favorite applica¬ 
tion, it calculates the radiative transfer in a contaminated ice 
slab overlaying an optically thick granular medium. The con¬ 
tamination in the slab can be of any type : ice, minerals or even 
bubbles. The matrix can be constituted as well of any ice. Thus 
our model can be applied on Earth with water ice, but also on 
Mars polar region covered with CO 2 ice [26], on icy bodies, such 
as Jupiter's moon Europa (water ice), Neptune's moon Triton 
or dwarf planet Pluto (N 2 ice). Other applications in biology or 
industry are possible, as soon as the optical constants of each 
material are known. 

We considered in all the calculations every wavelength in¬ 
dependently. Thus, a spectrum in any spectral range can be 
built by computing every wavelength contribution at very high 
spectral resolution. The final objective is the comparison of the 
simulation to actual data, for analysis purposes. This makes this 
approach suitable for any spectroscopic measurement of slabs 
(made of ice or other material), overlaying optically thick mate¬ 
rial (granular or other material), from laboratory to spatial probe 
measurement. For the planetary science case, these results will 
be down-sampled at the instrument's wavelength resolution, 
using its PSFs. 

One major hypothesis in this work is that we suppose an 
isotropic behavior of the inclusions. In the future, we plan to add 
the particles phase function to improve this point. We also plan 
to normalize the probability distribution function describing 
the roughness of the surface, to extend the applicability of the 
model. 
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NOTATIONS 

tx! phase angle (°) 

7 c compactness of the matrix : volume of matrix per 

unit of volume 

£ slope azimuth angle (°) 

£spec slope azimuth angle of a facet in specular condi¬ 
tions (°) 

9 roughness parameter (°) 

9 slope angle (°) 

9spec slope angle of a facet in specular conditions(°) 

©ip transmission factor of a type k inclusion 

© transmission factor of the slab containing inclusion 

under isotropic illumination 

©' transmission factor of the slab containing inclusion 

under collimated illumination 

v optical path (m) 

pp radius of a type k inclusion (m) 

crp geometrical cross section for type k inclusions (m 2 ) 

cr e p extinction cross section for type k inclusions (m 2 ) 

( cr e ) mean extinction cross section (m 2 ) 

cr s p scattering cross section for type k inclusions (m 2 ) 

(cr e ) mean scattering cross section (m 2 ) 

t optical depth of the matrix containing inclusions 

i p azimuth angle (°) 

co single scattering albedo of the matrix containing 

inclusions 

to s single scattering albedo of the granular substrate 

Qc Solid angle subtended by the sensor (sr) 

Qj Solid angle subtended by the source (sr) 

a m absorption coefficient of the matrix 

cijp absorption coefficient a type k inclusion 

a (9, £) probability of occurrence for the slope($, £) 

A Surface of a pixel (in 2 ) 

Af Surface of a facet (m 2 ) 

D thickness of the slab layer (m) 

D' apparent length of the first transit through the slab 

layer for one ray (m) 

D' mean apparent length of the first transit through 

the slab layer for one ray (m) 

e emergence angle (°) 

F incident power flux in the radiation's direction 

(W.m -2 ) 

i incidence angle (°) 

k m Imaginary part of the optical index of the matrix 


kjp Imaginary part of the optical index of a type k in¬ 
clusion 

L radiance (W.m -2 .sr -1 ) 

n m Real part of the optical index of the matrix 

rijp Real part of the optical index of a type k inclusion 

N number of facets within a pixel: N 1 

N S pec number of facets within a pixel satisfying specular 
conditions 

M total density of inclusions inside the matrix 
M k density of inclusions of type k inside the matrix 
p probability or probability per unit of length 

P power (W) 

Q sk scattering efficiency for type k inclusions 

r bidirectional reflectance (sr -1 ) 

Tjp diffusive reflectance for a type k inclusion 

r m diffusive reflectance of the matrix 

r s diffusive reflectance of the granular substrate 

tf Fresnel reflection coefficient r f — R 2 ± + R 2 

Rx Fresnel reflectivites for perpendicular polarization 

R|| Fresnel reflectivites for perpendicular and parallel 

polarization 

Ro reflection factor of the slab containing inclusion 
under isotropic illumination 

Rq reflection factor of the slab containing inclusion 

under collimated illumination 

/?' 7?" _ c' 

JXq s e 

R-Diff diffuse reflectance factor of a slab containing inclu¬ 
sion over a granular substrate 

R S pec specular reflectance factor of a slab containing in¬ 
clusion over a granular substrate 

Rf 0 f reflectance factor of a slab containing inclusion over 
a granular substrate 

S shadowing function 

S' e external reflection coefficient of the slab under colli¬ 
mated illumination 

S e external reflection coefficient of the slab under 
isotropic illumination 

S; internal reflection coefficient of the slab under 

isotropic illumination 

S ek external reflection coefficient of a type k inclusion 

Sjp internal reflection coefficient of a type k inclusion 

To transmission factor of the slab containing inclusion 

under isotropic illumination 

Tg transmission factor of the slab containing inclusion 

under collimated illumination 



